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Abstract 

We develop an algorithm for the computation of general Fourier integral operators associated with 
canonical graphs. The algorithm is based on dyadic parabolic decomposition using wave packets and 
enables the discrete approximate evaluation of the action of such operators on data in the presence of 
caustics. The procedure consists in the construction of a universal operator representation through the 
introduction of locally singularity-resolving diffeomorphisms, enabling the application of wave packet 
driven computation, and in the construction of the associated pseudo-differential joint-partition of unity on 
the canonical graphs. We apply the method to a parametrix of the wave equation in the vicinity of a cusp 
singularity. 

1 Introduction 

In this paper, we develop an algorithm for applying Fourier integral operators associated with canonical 
graphs using wave packets. To arrive at such an algorithm, we construct a universal oscillatory integral 
representation of the kernels of these Fourier integral operators by introducing singularity resolving diffeo- 
morphisms where caustics occur. The universal representation is of the form such that the algorithm based 
on the dyadic parabolic decomposition of phase space previously developed by the authors applies 0. We 
refer to 18] [TUJ [TT) for related computational methods aiming at the evaluation of the action of Fourier 
integral operators. 

The algorithm comprises a geometrical component, bringing the local representations in universal form, 
and a wave packet component which yields the application of the local operators. Here, we develop the 
geometrical component, which consists of the following steps. First we determine the location of caustics 
on the canonical relation of the Fourier integral operator. For each point on a caustic we determine the 
associated specific rank deficiency and construct an appropriate diffeomorphism, resolving the caustic in 
open neighborhoods of this point. We determine the (local) phase function of the composition of the Fourier 
integral operator and the inverse of the diffeomorphism in terms of universal coordinates and detect the 
largest set on which it is defined. We evaluate the preimage of this set on the canonical relation. We continue 
this procedure until the caustic is covered with overlapping sets, associated with diffeomorphisms for the 
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Figure 1: Projection A(y,rj2) of a slice £ = £o of the canonical relation A associated with a half- wave 
equation in the vicinity of a caustic (red solid line, the blue dashed lines indicate the neighborhood of the 
singularity) caused by a low velocity lens. The white solid lines are connected to a regular grid in x by bi- 
characteristics. The black dot indicates the center of an open neighborhood of conjugate points (#o, £o) ^ 
(2/0? Vo) f° r which the projection onto standard microlocal focal coordinates (y, £) is not diffeomorphic. 



corresponding rank deficiencies. Then we repeat the steps for each caustic and arrive at a collection of open 
sets covering the canonical relation. 

The complexity of the algorithm for general Fourier integral operators as compared to the non-caustic 
case arises from switching, in the sets covering a small neighborhood of the caustics, from a global to a local 
algorithm using a pseudodifferential partition of unity. 

As an application we present the computation of a parametrix of the wave equation in a heterogeneous, 
isotropic setting for long-time stepping in the presence of caustics. 

Curvelets, wave packets 

We briefly discuss the (co)frame of curvelets and wave packets l9l fT3ll24l . Let u £ L 2 (R n ) and consider its 
Fourier transform, u{£) = J u{x) exp[— £)] dx. 

One begins with an overlapping covering of the positive £1 axis (£' = £1) by boxes of the form 

n-l 

(1) 

where the centers as well as the side lengths L' k and L'£, satisfy the parabolic scaling condition 

e k ~2 k , L' k ~2 k , Ll~2 k '\ asfc^oo. 

Next, for each k > 1, let v vary over a set of ~ 2 k ^ n ~ 1 ^ 2 uniformly distributed unit vectors. Let denote 
a choice of rotation matrix which maps v to ei, and B v ^ = Q~\Bk- In the (co-)frame construction, one 

encounters two sequences of smooth functions on R n , \v,h an ^ $v,k, eacn supported in B u ^, so that they 
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form a co-partition of unity, Xo(OA)(0 + l>2k>i Xi/,fe(0^,fe(0 = and satisfy the estimates 

One then forms ^„,fc(f) = PjT 1/2 A/,fc(0> ^,fe(0 = P/T 1/2 X^,/e(0> witn Pfc = vol^), satisfying the 
estimates 

VAT: l ^' /c(x)l | <c N 2^ n+1 ^ A (2 k \{v,x)\^2 k ^\\x\\)- N . (2) 

To obtain a (co)frame, one introduces the integer lattice: Xj := (ji, . . . , j n ) G Z n , the dilation matrix = 

f - ( ^ ?), x r n_1 J , det L> fc = (27r)- n p fc , and points x v * = Q-^D^Xj. The frame elements 

n \ U n _ixi L k l n -i J 

(k > 1) are then defined in the Fourier domain as <£ 7 (£) = fiv,k(0 exp[— £)], 7 = ^, fc), and 

similarly for ^ 7 (£). The function is referred to as a wave packet. One obtains the transform pair 



J u(x)xpy(x) dx, u(x) = ' S y^^u 1 Lp 1 (x). 



(3) 



2 Fourier integral operators and caustics 

We consider Fourier integral operators, F, associated with canonical graphs. We allow the formation of 
caustics. 

2.1 Oscillatory integrals, local coordinates 

Let (y : xj i , £j. ) be local coordinates on the canonical relation, A say, of F, and Si a corresponding generating 
function: If, at a point on A, (dy, dxi) are linearly independent and dxj vanishes, then (dy, dxi, d^j) are 
coordinates on A nearby, /U J = {l,...,n}, I D J = {0}, and one can parameterize A as (Xj(y, xj, £j) — 
^j, where = Xj(y, xi, £j) locally on A (cf. (T8), Thm 21.2.18). The fact that a (possibly empty) 
set / exists follows from the canonical graph property, i.e. that (y, 77) are local coordinates and dy linearly 
independent. Then 

dSj dSj 

The coordinates are standardly defined on (overlapping) open sets Oi in A, that is, (y, xj. , £j. ) — » r(y, x/. , £j. ) 
is defined as a diffeomorphism on Of, let z = 1, . . . , N. The corresponding partition of unity is written as 

N 



5^A(r) = l, rG A. (5) 



2=1 

In local coordinates, we introduce 



riiy.xi.^j.) = r i (r(y,x Ii ,£ sJ .)). (6) 
Then (F<pJ(y) = E^ii^M with 

(Fi(p 7 )(y) = / / f^x/.^Ja^xj^OJ^ (7) 



The amplitude 0^(3/, x/. , ^ ) is complex and accounts for the KMAH index. 
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We let E^ denote the stationary point set (in 6) of <p = (j)(y, x,Q). The amplitude can be identified with a 
half-density on A. One defines the 2n-form on E^, 



d^ Ad 



d0 1 



A... Ad 



= dyi A... A dy n A dx\ A ... A dx n A d0\ A ... A d#Ar. 



In the above, we choose A = (y, x/, J^-) as local coordinates on A, while # = £j. Then we get 



d = |A |-VAiA...AdA 2 n|, 



dxjdxj d^jdxj 



-1; 



A is identified with (?/, x/, £j). The corresponding half-density equals | | — 1//2 |<iAi A ... A d\2 n \ 1 ^ 2 - 

Densities on a submanifold of the cotangent bundle are associated with the determinant bundle of the 
cotangent bundle. Let a® denote the leading order homogeneous part of a^. The principal symbol of the 
Fourier integral operator then defines a half-density, a^d l J 2 . That is, for a change of local coordinates, if the 
transformation rule for forms of maximal degree is the multiplication by a Jacobian j, then the transformation 
rule for a half-density is the multiplication by \j\ 1 ^ 2 . In our case, of canonical graphs, we can dispose of the 
description in terms of half-densities and restrict to zero-density amplitudes on A. 



2.2 Propagator 



The typical case of a Fourier integral operator associated with a canonical graph is the parametrix for an 
evolution equation fT4l[T5lL 



[d t + iP(t, x, D x )]u(t, x) = 0, u(t , x) = cp 7 (x) 



(8) 



on a domain X C W 1 and a time interval [to, T], where P(t, x, D x ) is a pseudodifferential operator with 
symbol in S\ ; we let p denote the principal symbol of P. 

For every (x, £) G T*X\{0}, the integral curves (y(x, £; t, to), 77(2;, £; t, to)) of 



dp(t,y,rj) (h]_ 

dt 



dt 



dp(t,y,rj) 
dy 



(9) 



with initial conditions to, to) = x and 77(x, to, to) = £ define the transformation, from £) 

to (2/, 77), which generates the canonical relation of the parameterix of ([5]), for a given time t = T; that is, 

The perturbations of (y, rj) with respect to initial conditions (x, £) are collected in a propagator matrix, 



nOr,£;t,t ) 



Wi w 2 
W 3 W A 



d x y d^y 
d x rj d^r] 



(10) 



which is the solution to the 2n x 2n system of differential equations 



dU 

~dt 



( d 2 p 



dr]dy 
d 2 p 



d 2 p 

dr]dr] 
d 2 p 



\ 



(t,y,v) 



n(x,^;t,t ), 



(11) 



\ dydy y ' u dydr] v 
known as the Hamilton- Jacobi equations, supplemented with the initial conditions l25ll26l 



n(x,^;t ,t ) 



I 
I 



(12) 
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Figure 2: Illustration of canonical relations x ( to P) an d X (bottom) of operators F and F associated with 
a half- wave equation: (bi-)characteristics ("rays") in y for initial conditions (x2 = #2,0? £ = £0) and 
(£2 = ^2,0 7 £ = £0), respectively, for evolution through a low velocity lens (see Section [5]). The black 
circles on the left indicate the conjugate points corresponding to the initial conditions. 

Away from caustics the generating function of A is S = S(y, £) (Ii — 0), which satisfies 



d2s ( t\ 
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= Wf 
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(13) 
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(15) 



upon substituting x = x(y,£]to,T)) denoting the backward solution to ^ with initial time T, evaluated at 
to. The leading-order amplitude follows to be 

a(y, C/ICI) = y/V det e/ICh *o, T), T, i ), (16) 

reflecting that a is homogeneous of degree in £. 

In the vicinity of caustics, we need to choose different coordinates. Admissible coordinates are directly 
related to the possible rank deficiency of W\ : One determines the null space of the matrix W\ and rotates 
the coordinates such the null space is spanned by the columns indexed by the set 1{. Then (y, xj i , form 
local coordinates on the canonical relation A, as in the previous subsection, and Oi is given by the set for 
which the columns indexed by Ii span the null space of W± . 

3 Singularity resolving diffeomorphisms 

We consider the matrix Wi(x(y, £; to, T), £; T, to) for given (T, to) at yo = 2/(#o> £o; ^, to) and £ = £o and 
determine its rank. Suppose it does not have full rank at this point. We construct a diffeomorphism which 
removes this rank deficiency in a neighborhood of ro = (yo-, Vo'i £o) £ A, where 770 = 77 (#0? £0; to)- 

To be specific, we rotate coordinates, such that £0 = (1? 0, . . . , 0) (upon normalization). Let us assume 
that the row associated with the coordinate X2 generates the rank deficiency. (There could be more than one 



row / coordinate.) We then introduce the diffeomorphism, 

Q:x^x={x 1 - ^(x 2 - Oo)2) 2 ,^2, • • • ,x n )\ 
to preserve the symplectic form, we map 

(6,6 + ^(^2 - (^0)2) ^1, 6, • ..,fn), 

yielding a canonical transformation Cq : £) (#,£). We note that Cq (2:0 , £0) = (^o?^o)- The 
diffeomorphism Q can be written in the form of an invertible Fourier integral operator with unit amplitude 
and canonical relation given as the graph of Cq (see Appendix |A|. 

The canonical transformation, Cq 1 , associated with Q~ x is given by 



x —> x — (xi 



2 {X2 



xoa) 2 ,x 2 , ■ ■ .,x n ), 



We introduce the pull back, Q*u(x) = u(Q 
ing propagator matrices are given by 

/ 1 






(6,6- 



a(x 2 - 
= u(xi - 



^0,2)fl, • • • ,in)- 




-a(x 2 - 2:0,2) 
1 




dx dx 



n; 



dx 








a(x 2 - 2:0,2) 
1 









V 



(X ) 2 ) 2 ,^2,... 



1 

a(x 2 — 2:0,2) 








1 

-a(x 2 — 2:0,2) 




x n ). The correspond- 



(17) 



(18) 
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which are easily verified to be symplectic matrices. In the more general case, each coordinate xj generating 



a rank deficiency yields additional non-zero entry pairs l^ 1 , H 1 , and ^ in the above 



dX! ggj ggj 



0& 



9fi 



propagator matrices. It follows that the composition (x, £ ) ^ (2:, £) ^ (y, 7/) generates the graph of a 
canonical transformation, x say, which can be parametrized by (y, £) locally on an open neighborhood of 
(2/0, £(#0? £o))- We denote the corresponding generating function by S = S(y, £). We can compose F with 
Q _1 as Fourier integral operators: F = FQ~ X . The canonical relation of F is the graph of x- m summary: 
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Figure 3: Caustic surfaces E(y, £) (dark gray) and E(y : £) (light gray) of A and A corresponding to propa- 
gation through a low velocity lens (cf. Section [5}: The singular regions of A and A do not intersect. 

For each given type of rank deficiency (here, in x 2 ) and each (xq, £0) within this class, there is an open set 
O(x ,£o) on which the coordinates (J, J) are valid. These sets form an open cover, and we obtain a family of 
diffeomorphisms parametrized by (xq, £0); there exists a locally finite subcover, and we only need a discrete 
set to resolve the rank deficiencies everywhere. We index these by j = 1, . . . , N{ and construct a set of 
diffeomorphisms, {Qij}f=i, which resolve locally the rank deficiency leading to coordinates (y, xi i , 
We write 

(y.x^^j,) ^ (y,i) 

c Qi . 

A 3 r = (y,r]\x,£) -A (j/,77; x, f) = f e A^- 

d 2 S- 

We write 0{ for the image of 0{ under the diffeomorphism on the level of Lagrangians. Let the matrix 



dyd£ 

in the above be nonsingular on the open set Uij, and introduce 0^ = Uij fl 0$ C A^ . This set corresponds 
with a set Oij C A. We subpartition Oi = U J= i v .. ) Ar i Oij. The corresponding partition of unity now reads 

N Ni 

X)E^W = while ^(2/,^, = A;(r(y, j = l,...,JVi. (19) 

i=i j=i 

Then (F<pJ(y) = ^=1 EfiiWi^)^) with 



( F ij¥i)(y) = J J Aj(y,^,OJ«i(y^/ i ,OJexp[i(S' i (?/,x / .,OJ - (Oi^Ji)] ^7( x ) dxd Oi- (20) 
Inserting the diffeomorphisms, we obtain 

(F ij( p y )(y) = J Aijiy,!) eMiSij(y,£)} Qj^(l) d£ (21) 
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Figure 4: Iso-amplitude surface of the partition functions F$ i = 1,3, associated with Q i = 1 

(left): the joint admissible set 0\ U O3 comprises the exterior of the two sheets. Iso-amplitude surface of 
Tij(x(x, £), £)) for £ = = and a = 1 (z = 2. j = 1) (right): the admissible set 

contains the region on the back of the sheet. Bottom: boundaries dOi, i = 1, 3 (dashed curves) and dOij, 
i = 2, j = 1 (solid curves) of the admissible domains: Clearly, the joint admissible set 0\ U O3 U O21 covers 
A. 

The amplitude (y, £) and phase function Sij (y, £) — (£, x) are obtained by composing Fij with Q" 1 
as Fourier integral operators and changing phase variables. It is possible to treat this composition from a 
semi-group point of view. Then, to leading order, we get 



where 
in which 



Aij(y,i) = fij(y,i) cLij(y,i), 
Ai(y,0 = Ai(r(y,0), 

Ai(r(r)) = Ai(r). 



(22) 

(23) 
(24) 



Moreover, a^- (y, £) can be obtained as follows. If II is the propagator matrix of the perturbations of x, then 
the propagator matrix of the perturbations of x is given by: = II IIq* . Then 



\ 



l/det 



dSij(y,i) 



(25) 
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Tij(x(x, £), £) = £o) (center), and the partition function Tij(x(x, £),£(x, £) = 1.67) for Oij realizing 
the partition of unity with Ti(x, £ = £o)- 



where det J * s obtained as the determinant of the upper-left sub-block of fl^ . To accommodate 

a common notation, we set Qij = I (Ni = 1) if Ii = 0. In the further analysis, we omit the subscripts ij 
where appropriate. 



Expansion of the cutoff functions 

The application of our algorithm involves the re-decomposition of Q*(p 7 into wave packets. The key novelty 
is constructing a separated representation of the partition functions. 

Consider our oscillatory integral in (y, £) including the cutoff r(y, £). P(y, £) is homogeneous of degree 
zero in £ and is a classical smooth symbol (of order 0). We "subdivide" the integration over £. A possible 
procedure involves obtaining a (low-rank) separated representation of P(y, £) on the support of each relevant 
box in £013111, 

A2/,0 = £f?(2/)f^(0, feB^. (26) 

/3=1 

(Basically, this can be obtained using spherical harmonics in view of the fact that the £ is implicitly limited 
to an annulus.) One can view this also as windowing the directions of <f into subsets (cones) using f 2 (0 an d 
then constructing f ^ (y) according to the smallest admissible set in y for the /3-range of directions. 
The oscillatory integral becomes 

Jv,k „ 

(F^)(y) = J2J2 f i^ / a(y,v)exp[iS(y,i)}Tl(i) |x^(£)| 2 Q*^(0df. (27) 

/3=1 ^ 

One can view (£)Xv,k(£) as a subdivision of the box B v ^. We know that | J v ^\ 1 as fc 00 since the 
cone of directions in B v ^ shrinks as y/k. Hence, for large k this does not involve any action. 

The procedure allows a subdivision for coarse scales, as long as the scaling is not affected for large k. If 
the subdivision is too "coarse" then parts of the integration will be lost. 



4 Computation 

Here, we develop an algorithm for applying Fourier integral operators in the above constructed universal 
oscillatory integral representation. The algorithm makes use of the wave-packet based "box- algorithm" com- 
putation of the action of Fourier integral operators associated with canonical graphs in microlocal standard 
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Figure 9: Illustration of diffeomorphism Q, Q~ x and re-decomposition for a wave packet (p 7 (x) at frequency 
scale k = 2. Top row: (f 7 (x) (left) and pull-back Q*p 7 (x) (right). Rows 2 to 4: Re-decompositions 
k Uj(fj(x) of Q*ip 1 (x) using 3, 7, and 9 boxes B„^, respectively (right column), and the corresponding 

image (^2$ k u^(p^j (x) under the action of Q _1 (left column). An insufficient number of boxes 

alters the amplitudes and the minimum phase property of wave packets. Increasing the number of boxes 
yields satisfactory results in an open neighborhood of (#o, £o)- 
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focal coordinates (y, £) 13. It is based on the discretization and approximation, to accuracy 0(2 fc / 2 ), of 
the action of on a wave packet <£j, 



(^^)(i/)«^(y^)E a S(y) E ei 



x^(£)l 2 4;l(£) 



(28) 



The procedure relies on truncated Taylor series expansions of Sij(y,£) and A(y,£) near the microlocal 
support of tpj, along the v = ^'/l^'l ax * s anc * m me £" directions perpendicular to the radial £ = £' direction. 
Here, T„^{y) is the backwards-solution 



_ dSij(y,u) 



and a~ r ^,(y) and #~ r j|.(0 are functions realizing, on a separated tensor-product representation of the 

slowly oscillating kernel appearing in the second-order expansion term of Sij , 



exp 



2 f/ \ d^ 2 



(29) 



We construct the functions o^jj, (?/) and $ ~ r j|, (£) and the tensor product separated representation from prolate 
spheroidal wave functions |[6j EZl E3 [23 • We refer to El for a detailed description of the box- algorithm. 



Preparation step 

We begin with determining the sets Oi by computing the perturbations of the integral curves (y(x,£),r](x,£)) 

dy_ 

dx ' 



with respect to initial conditions (x, £) and monitoring the null space of the matrix 1^, as detailed in Section 



2.2 



For parametrices of evolution equations, this involves evaluation of the propagator matrices n(x, £). 
Then, for each set C^, we detect Uij (and consequently O^j) in a similar way, as the sets on which the upper 
left sub-block of fl^ = II Ilg. . has full rank. Here Ilg 1 (x, £) is given by ([18]). 

We then proceed with the construction of the partition of unity. Since the partition functions enter the 
computation as pseudodifferential cutoffs in the construction of the amplitude (cf. ([22])), requiring the back- 



wards solutions x(y, £) (compare ( 13 -16)), we perform our numerical construction in coordinates (x, £). We 
obtain Tij(y, £) upon substituting ?/ = 2/(5, implied by the canonical relation Xij- F° r the construction of 
the partition functions , we choose double-exponential cutoffs of the form 

exp(— exp(d(x, £)) 

mimicking a Cq° cutoff, with appropriate normalization and truncated to precision e. Here d(x, £) is a 
function measuring the distance of the point from the boundary d?7^ of the set Uij(x,£). The 

partition of unity is then formed by weighting f^-(x, £) on the overlaps of the sets Uij(x,£) such that 
Y^ij f ij(f(x, £)) = 1. Finally, we construct the separated representations in coordinates by win- 

dowing the directions of £ into subsets using f % (£), realizing a subdivision into £ cones. 

Diffeomorphism 

We evaluate Q in the Fourier domain. The data </? 7 (x) enter the box algorithm via the coefficients u 7 of their 
discrete almost symmetric wave packet transform fT3lL allowing the fast evaluation of the Fourier transform 
of the data at a set of frequency points <^' fc limited to the box B v ^. From these, we obtain Q*-ip 1 {x) via 
evaluation of adjoint unequally spaced FFT CUE) at points x(x). 



14 



Application of the box algorithm 

We are now ready to compute the action {F i jip 1 ){y) (cf. ^2A\ ) by applying the box algorithm (cf. ( [28] )) to the 
pull-back Q*j(f 7 (x). First, we compute the discrete almost symmetric wave packet transform of Q*jp 7 (x), 
yielding its wave packet coefficients Note that numerically significant coefficients are contained 

in a small set of boxes neighboring the direction £ = £o/|£o|- We subdivide each box according to the 
separated representation of f ^ (cf. ([27])). Then, we apply the box algorithm to each subdivision, indexed 
by triples (/?, z>, k), f3 = 1, . . . , J^. Here, the Taylor series expansion of the generating function Sij(y, £) 
underlying the box algorithm is constructed about the central £ direction within the support of (OXf?,fe (0» 
accounting for the induced subdivision of the box By^. Note that sub-dividing into | cones results in a 
reduction of the range of £ orientations in each element z>, fc) of the subdivision, as compared to the £ 
range contained in This reduces the number R of expansion terms necessary in the separated tensor 

product representation for yielding prescribed accuracy, and effectively counter-balances the increase by a 
factor Jf^fc, evoked by the separated representation of T^ , of the number of times the box-algorithm has to 
be applied. 

Operator hierarchy 

The operators F^ for which Qij = I, say, are directly associated with the canonical relation and 
involve only computations on Ap. In the algorithm, we reflect this physical hierarchy of the operators F^ 
in the construction of the partition of unity. First, we construct a partition of unity for these hierarchically 
higher operators. Then, we construct a joint partition of the remaining operators on the sets which are not 
covered by the sets for which = L 

Re-decomposition 

Starting from a single box B u ^, re-decomposition of Q*tp 1 (x) results in a set of boxes B„^ yielding nu- 
merically non-zero contribution to the solution. The number of boxes entering the computation is directly 
proportional to the computational cost of the algorithm. In applications, we therefore aim at keeping this 
number small and consider only a subset of boxes, yielding the most significant contributions. We choose 
this subset such that on an open neighborhood of (#o, £o) 

QTj-Qa « i- 

We can estimate the energy loss induced by the restriction to subsets of and re-normalize the solu- 
tion. We illustrate the impact of choices of subsets containing different numbers of boxes on the numerical 
accuracy of the diffeomorphic identity in Fig. [9] 

Furthermore, the re-decomposition of Q*(p 7 (x) yields in general, under the action of ^-values 
outside the set B u ^, £(#,£) D B v ^. We monitor and do not consider their contribution in our 

computation if \x u ,k(£(%i 0)1 ^ s below a given threshold. 

5 Numerical example 

We numerically illustrate our algorithm for the evaluation of the action of Fourier integral operators associ- 
ated with evolution equations. We consider wave evolution under the half- wave equation, that is, the initial 
value problem ^ with symbol 

in n = 2 dimensions. Here c(x) stands for the medium velocity. 
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Figure 10: Illustration of operator action on a wave packet (p 7 (x) at frequency scale k = 2: Contribution 
of operators Fi (i = 1, 3) associated with = I (top left), contribution of operator (i = 2, j = 1) with 
diffeomorphism parameters (£o = tt/2, £2,0 = 7 ^ = 1)> resolving the singularity in the tip of the caustic 
(top right), and joint action of Fi and (bottom left). Time domain finite difference reference (bottom 
right). In the operator computation, we consider 9 boxes and a separated representation with J v ^ = 1 
term. 



Heterogeneous, isotropic model 

We choose a heterogeneous velocity model 

c(x) = c + ftexp(-|x - x | 2 /cr 2 ), 

containing a low velocity lens, with parameters Co = 2km/ s, k = —OAkm/s, a = 3km, and xq = 
(0, 14) km. As the initial data, we choose horizontal wave packets at frequency scale k — 2 and k = 3, 
respectively, in the vicinity of the point x' = (0, 5) km. We fix the evolution time to T = 7s. With this choice 
of parameters, most of the energy of the solution is concentrated near a cusp-type caustic. We illustrate the 
induced sets Oi and the joint partition of unity I\ in Fig. [4]and[5] 



Operator factorization 

We partition the Lagrangian A into three sets 0^, i = {1,2,3}. The sets i = {1,3} are separated by 
the caustic. For these sets, we can choose coordinates (?/,£)> nence Q% — I- The set i = 2 contains the 
caustic. For illustration purposes, in the factorization Fij of Fi for i = 2, we choose to compute the operator 
j = 1, which resolves the singularity in an open neighborhood of the point indicated by a black dot on the 
Lagrangian plotted in Fig. [T] This neighborhood contains the cusp of the caustic. Furthermore, we limit our 
separated representation to one term, J u k = 1 (for the corresponding admissible sets and partition functions, 
see Fig. |6||8](top rows)). We restrict the computation of F^ for the initial data at frequency scale k = 2 
(k = 3) to 9 (11) boxes neighboring the v direction, respectively. 



Results 

In Fig. [10| we plot the contributions of the different components in the factorization of the propagator acting 
on a single horizontal wave packet at frequency scale k — 2, and compare to a time domain finite difference 
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Figure 11: Left column: Joint contribution of the operators Fi and F^ acting on a wave packets ip 7 (x) at 



frequency scale k — 3 (compare Fig. 10 (bottom left) for a wave packet at frequency scale k = 2). Center 
column: Time domain finite difference reference. Right column: The equivalent to the left column when 
using a separated representation with J v ^ = 11 terms (note that for computational reasons, only 1 box 
has been used in the numerical evaluation of Fi and F^ with J v ^ = 11 terms). 
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computation. The support of the wave packet within the joint admissible set of the chosen factorization is 
mostly covered by the set O^, such that most of its energy is contributed by the operator F^ , for which 

Qij + I 

We observe that in the joint admissible set, our algorithm has effectively removed the singularity. We note 
that the phase of the operator computation matches the phase of the finite difference reference. This includes 
the KMAH index, which is best observed for operator F3, which exclusively contributes to the region beyond 
the caustic (cf. Fig. [10| top left). Furthermore, note that the amplitude obtained by our algorithm is slightly 
weaker than the true amplitude. This is consistent with the observations and discussion following Fig. [9] 
and results from the energy leakage induced by restricting the number of boxes in the re-decomposition step 
following the application of Q. We can compensate and re-normalize the amplitude by monitoring the energy 
loss resulting from the restriction (in Fig. [10| we have not re-normalized the amplitudes). Finally, we note 
that our algorithm yields the correct result in an open neighborhood in the vicinity of the tip of the caustic, for 
which we have designed the operator Fij . In consistency with this fact, it is ineffective for yielding the image 
of the entire wave packet which, at this low frequency scale, has support extending beyond the admissible 
set of the operator factors we compute. 



These observations are further illustrated in Fig. 1 1 (left column), where we plot the contributions of 
the different components in the factorization of the propagator acting on horizontal wave packets, at higher 
frequency scale k = 3, centered at locations in the vicinity of the caustic tip (results of a time domain 
finite difference reference computation are plotted in Fig. [TT] (center column)). With these initial data, we 
explore the open neighborhood about the point for which the operator composition with resolves the 
singularity. Indeed, at this frequency scale, we can obtain the image of an entire wave packet with only one 
operator factor (cf. Fig. 11 (second row)). For the wave packet located slightly further above the tip of 
the caustic (top row), we observe a phase artifact in the region of overlap of Oi=\ and Oij, which can be 
explained as follows: The restriction of the separated representation for Fi to one term only induces that the 
computation of the geometry (bi-characteristics) for the entire box B v ^ is exclusively based on one single 
direction v. This results in inaccuracies in regions close to the caustics where slight perturbations in £ yield 
large variations in y. Furthermore, as discussed above, wave packets exploring the regions beyond the tip of 
the caustics eventually start to leave the admissible set for F^ (third and bottom line). 

We note that both for removing the phase artifact of Fi close to the caustic, and for enlarging the admissi- 
ble set, it is necessary to increase the number of terms J v ^ in the separated representation ( [26] ) (compare Fig. 



[8]). This is illustrated in Fig. 11 (right column), where the joint contributions of operators Fi and F^ with 
J Vj k = 11 terms in the separated representation are plotted. Here, the expansion functions f f (0 are con- 
structed as as cones in £ with a squared cosine cutoff window. For practical reasons and illustration purpose, 
only one single frequency box has been used in the computation (cf. Fig. [9]). While this restriction to 
only one frequency box affects the amplitudes and the phases in the tails of the wave packet, the separated 
representation remains nonetheless effective in resolving the issues observed above: the admissible set is 
extended beyond the caustic and the inaccuracies in the regions of overlap of sets Oi and as well as in 
the regions close to the caustics are considerably reduced. 



6 Discussion 

We developed an algorithm for the evaluation of the action of Fourier integral operators through their fac- 
torization into operators with a universal oscillatory integral representation, enabled by the construction of 
appropriately chosen diffeomorphisms. The algorithm comprises a preparatory geometrical step in which 
open sets are detected on the canonical relation for which specific focal coordinates are admissible. This 
covering with open sets induces a pseudodifferential partition of unity. Then, for each term of this partition, 
we apply a factorization of the associated operators using diffeomorphisms reflecting the rank deficiency 
and resolving the singularity in the set. This factorization admits a parametrization of the canonical graph 
in universal (y, £) coordinate pairs and enables the application of our previously developed box algorithm, 
following the dyadic parabolic decomposition of phase space, for numerical computations. Hence, our algo- 
rithm enables the discrete wave packet based computation of the action of Fourier integral operators globally, 
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including in the vicinity of caustics. This wave packet description is valid on the entire canonical relation. It 
can now enter procedures aiming at the iterative refinement of approximate solutions, and drive the construc- 
tion of weak solutions via Volterra kernels fTHT2ll. 

In the special case of Fourier integral operators corresponding to parametrices of evolution equations, for 
isotropic media, an alternative approach for obtaining solutions in the vicinity of caustics has been proposed 
previously OQjOEO). It consist in a re-decomposition strategy following a multi-product representation of 
the propagator. Here, we avoid the re-decompositions and operator compositions following the discretiza- 
tion of the evolution parameter, reminiscent of a stepping procedure. What is more, our construction is not 
restricted to parametrices of evolution equations, but is valid for the general class of Fourier integral op- 
erators associated with canonical graphs, allowing for anisotropy. The cost of the algorithm resides in the 
construction and application of the separated representation of the pseudodifferential partition of unity. 



A Fourier integral representation of Q and Q 1 

We write (Q*u)(x) = u(X(x)), ((Q- 1 )*£t)(x) = u(X(x)). That is, X = Q' 1 and X = Q. The 
diffeomorphisms Q and Q~ x define the Fourier integral operators with oscillatory integral kernels, 



A Q (x,x) = J e-^-^^dC, A Q -i(x,x) = J G -Ht>2-*(* 
The generating functions are 



(30) 



S Q (x,0 = (£,X(x)), S Q -i(x,£) = (lX(x)), 
respectively. The canonical relations are the graphs of Cq and Cq-i, and are given by 
A Q = {(x = X- 1 (x),{Z,d £ X)\ ii=x -i {x y,x,Z)}, A Q -i ={(x = X-\x),{id x X)\ x= x- 1(£ y,x,i)}. 
The Hessians yield a unit amplitude: 



det 



d 2 (^X(x)) 



dxd^ 



det 



d 2 (lX(x)) 



dxd^ 



1. 



Substituting the particular diffeomorphism, we obtain: 



=X-i(5) 



{ld x X)\ x=jt - H , 



( 1 -a(x 2 - x , 2 ) 
1 

1 

V ! ! : 

| 2 - a(x 2 - x ,2)ii 



References 

[1] F. Andersson, M. V. de Hoop, H.F. Smith, and G. Uhlmann. A multi-scale approach to hyperbolic 
evolution equations with limited smoothness. Comm. Partial Differential Equations, 33:988-1017, 
2008. 

[2] F. Andersson, M.V. de Hoop, and H. Wendt. Multi-scale discrete approximation of Fourier integral 
operators. SIAM Multiscale Model Simul, 10:111-145, 2012. 



19 



[3] G. Bao and W.W. Symes. Computation of pseudodifferential operators. SI AM J. Sci. Compute 17:416- 
429, 1996. 

[4] G. Beylkin, V. Cheruvu, and F. Perez. Fast adaptive algorithms in the non-standard form for multidi- 
mensional problems. Appl. Comput. Harmon. Anal, 24:354-377, 2008. 

[5] G. Beylkin and M.J. Mohlenkamp. Algorithms for numerical analysis in high dimensions. SI AM J. Sci. 
Comput, 26(6):2133-2159, 2005. 

[6] G. Beylkin and K. Sandberg. Wave propagation using bases for bandlimited functions. Wave Motion, 
41:263-291,2005. 

[7] B. Bradie, R. Coifman, and A.Grossmann. Fast numerical computations of oscillatory integrals related 
to accoustic scattering, i. Appl. Comput. Harmon. Anal, l(l):94-99, 1993. 

[8] E. Candes and L. Demanet. Curvelets and Fourier integral operators. C. R. Acad. Sci. Paris, I(336):395- 
398, 2003. 

[9] E. Candes, L. Demanet, D. Donoho, and L. Ying. Fast discrete curvelet transforms. SIAM Multiscale 
Model. Simul, 5(3):861-899, 2006. 

[10] E. Candes, L. Demanet, and L. Ying. Fast computation of Fourier integral operators. SIAM J. Sci. 
Comput, 29(6):2464-2493, 2007. 

[11] E. Candes, L. Demanet, and L. Ying. A fast butterfly algorithm for the computation of Fourier integral 
operators. SIAM Multiscale Model. Simul, 7:1727-1750, 2009. 

[12] M.V. de Hoop, S.F. Holman, H.F. Smith, and G. Uhlmann. Regularity and multi-scale discretization 
of the solution construction of hyperbolic evolution equations with limited smoothness. Appl. Comput. 
Harmon. Anal, 33:330-353, 2012. 

[13] A. Duchkov, F. Andersson, and M.V. de Hoop. Discrete almost symmetric wave packets and multi-scale 
representation of (seismic) waves. IEEE T. Geosci. Remote Sensing, 48(9):3408-3423, 2010. 

[14] A. Duchkov and M.V. de Hoop. Extended isochron rays in shot-geophone (map) migration. Geophysics, 
75(4): 139-150, 2010. 

[15] A. Duchkov, M.V. de Hoop, and A. Sa Baretto. Evolution-equation approach to seismic image, and 
data, continuation. Wave Motion, 45:952-969, 2008. 

[16] A. Dutt and V. Rokhlin. Fast Fourier transforms for nonequispaced data. SIAM J. Sci. Comput., 
14(6): 1368-1393, 1993. 

[17] A. Dutt and V. Rokhlin. Fast Fourier transforms for nonequispaced data II. Appl. Comput. Harmon. 
Anal, 2:85-100, 1995. 

[18] L. Hormander. The Analysis of Linear partial Differenial Operators, volume IV. Springer- Verlag, 
Berlin, 1985. 

[19] H. Kumano-go, K. Taniguchi, and Y. Tozaki. Multi-products of phase functions for Fourier integral 
operators with an applications. Comm. Partial Differential Equations, 3(4):349-380, 1978. 

[20] J.H. Le Rousseau. Fourier-integral-operator approximation of solutions to first-order hyperbolic pseu- 
dodifferential equations i: convergence in Sobolev spaces. Comm. PDE, 31:867-906, 2006. 

[21] D. Slepian. Prolate spheroidal wave functions, Fourier analysis and uncertainty-IV: extensions to many 
dimensions, generalized prolate spheroidal wave functions. Bell Syst. Tech. J., November: 3009-3057, 
1964. 



20 



[22] D. Slepian. On the symmetrized Kronecker power of a matrix and extensions of Mehler's formula for 
Hermite polynomials. SI AM J. Math. Anal, 3:606-616, 1972. 

[23] D. Slepian. Prolate spheroidal wave functions, Fourier analysis and uncertainty-V: the discrete case. 
Bell Syst. Tech. 57:1371-1430, 1978. 

[24] H. Smith. A parametrix construction for wave equations with C 1,1 coefficients. Ann. Inst. Fourier, 
Grenoble, 48:797-835, 1998. 

[25] V. Cerveny. Seismic ray theory. Cambridge University Press, Cambridge, UK, 2001. 

[26] V. Cerveny. A note on dynamic ray tracing in ray-centered coordinates in anisotropic inhomogeneous 
media. Stud. Geophys. Geod, 51:411-422, 2007. 

[27] H. Xiao, V. Rokhlin, and N. Yarvin. Prolate spheroidal wave functions, quadrature and interpolation. 
Inverse problems, 17:805-838, 2001. 



21 



